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ABSTRACT 

High-redshift quasars at z > 6 have masses up to ^ 10® M 0 . One of the pathways to 
their formation includes direct collapse of gas, forming a supermassive star, precursor 
of the black hole seed. The conditions for direct collapse are more easily achievable 
in metal-free haloes, where atomic hydrogen cooling operates and molecular hydrogen 
(H 2 ) formation is inhibited by a strong external UV flux. Above a certain value of 
UV flux (Jcrit), the gas in a halo collapses isothermally at 10 '^ K and provides the 
conditions for supermassive star formation. However, H 2 can self-shield, reducing the 
effect of photodissociation. So far, most numerical studies used the local Jeans length 
to calculate the column densities for self-shielding. We implement an improved method 
for the determination of column densities in 3D simulations and analyse its effect on 
the value of Jcrit- This new method captures the gas geometry and velocity field 
and enables us to properly determine the direction-dependent self-shielding factor of 
H 2 against photodissociating radiation. We find a value of Jcrit that is a factor of two 
smaller than with the Jeans approach 2000 J 21 vs. ~ 4000 J 21 ). The main reason 
for this difference is the strong directional dependence of the H 2 column density. With 
this lower value of Jcrit, the number of haloes exposed to a flux > Jcrit is larger by 
more than an order of magnitude compared to previous studies. This may translate 
into a similar enhancement in the predicted number density of black hole seeds. 

Key words: black hole physics - methods: numerical - galaxies: formation - early 
Universe 


1 INTRODUCTION 


Observations of quasars at high redshifts indicate that su- 
permassive black holes (SMBHs) of several billion solar 
masses were already assembled in t h e firs t billion years 


after the big bang I 

Fan et al.l l2003l. 20061: Willott et al. 

201c 

: IVenemans et al. 

2OI3I: De Rosa et al. 2014: IWu et al. 

2015 

). The current record holders are a bright quasar, which 


hosts a SMBH with a mass of 2 x lO^M© at z = 7.085 
dMortlock et al.ll201lll . corresponding to ~ 800 million years 
after the b ig bang, a nd a SMBH with 1.2 x 10^°M© at 
2 = 6.30 JWu et al.l I2OI5II . It is still unclear how these 
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objects were able to acquire so much mass in this short 
period of time, which in turn raises questions about the 
formation mechanism and the involved physical processes. 
One pos sible explanation is the collapse of dense stellar 
clusters (|Portegijes_Z wartetaL _ 20041: lOmukai etTI] I 2 OO 8 I : 


iDevecchi &: Volonterill2009l : iLupi et al.ll2014ll . Another sce¬ 
nario involves stellar mass seed black holes with masses 
up to a few hundred M© that are the remnants of Pop¬ 
ulation III (Pop II I) stars and then grow by mass ac¬ 
cretion or mergers dMadau fc ReesI 2001 IIaiman_&Loed 


200l|_ VMonterie^_ah 2003 Yoo^jMiraldacEscud^200^ 


H.ajm ^~ 200j^ Pelui)essy_et alJ ' 2007: Tanaka fc Haim^ 


2OO9I : Whalen fc Frverll2012l : Madau et al.l 2014 ). However, 
already a simple order of magnitude argument shows that 
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2 T. Hartwig et al. 


this process involves some difficulties. Assuming accretion 
at th e Eddington limit, the e -folding time is 50 million 
years dMilosavlievic et al.ll2009l l. In 800 million years, a seed 
black hole accreting at the Eddington limit can therefore 
grow by a factor of ~ 9 x 10®, and so to reach 

a mass of 2 x 10® M© hy z = 7.085, it is necessary to 
start with a seed mass of ~ 200 M©. This is significantly 
larger than t he ma ss of a typical Pop I II stellar remnant 
dClark et ahl I2011alf9: iG reif et al ] l201ll: fs tacv et al.l l2012l : 


iLatif et al.l 2013al : Hirano et al. 20141 : Hartwig et al.ll2015ll . 

but not yet completely ruled out. Moreover, it is still an open 
question, how these high gas accretion rates could be sus¬ 
tained during the g rowth of the SMBH ( Alvureze^jjJ 200^ 
Mlffisavlievid et ^ l2009l : Ijohnson et al.l '20131 : ^on et al] 


2014h . 


A more promising formation scenario is the direct col- 


seed masses of — 

10®M,7i llReesI 19841: Loeb & Rasid 

1994 

:lBromm Loebll2003l:lBee:elman et alJl2006l:IVolonteril 

201 c 

: Shang et alJ 201C 

: Schleicher et al. 

20101: IChoi et al.1 

201 a 

:lLatif et al.1 

2013bl: 

Regan et al. 

2014 

Latif et al. 20141: 

Sugimura et al. 

2014 

Visbal et al. 


20L 

1: Agarwal et al.1 

2 OI 4 I: Latif et al. 

20151: 

Becerra et al. 

2 OI 5 I). To form such 


a massive se ed, a high mass inflow rate of > 0.1 M© yr 
equired dBegelman 2010l: Hosokaw a et ej'l l2012l . [201^ : 


ISchleicher et al.f 20131 : Ferrara et al. 20141 '). Sufficient condi¬ 
tions for such high mass inflow rates are provided in haloes 
with Tvir > 10“* K in which gas fragmentatio n and star 
format ion are suppressed during the collapse dLatif et al.l 
l2013bl ). To avoid fragmentation, the gas has to be metal 
free and a strong radiation background has to photodisso- 
ciate molecular hydrogen, which otherwise acts as a strong 
coolant. Under these specific conditions, the gas can only 
cool by atomic hydrogen and collapses monolithically to 
form a super massive star (SMS), which later on forms a 
SMBH seed dBegelmard l20ld: iHosokawa. et aid l2012l . l2013l : 
llnavoshi et al.ll2014l : Inavoshi fc Haiman 2014h or a quasi¬ 
star, which forms a stellar mass black hole that grows 


by swallowing its envelope dB egelrnan et al.l l2006l . 120081 : 
iBall et ^l201ll : ISchIeicher et al.l 2013ll . We will refer to this 
specific type of direct collapse as ‘direct collapse scenario’ 
hereafter. 

Based on the strength of the photodissociating radi¬ 
ation, the cloud either monolithically collapses close to 
isothermality, or is able to efficiently cool and to fragment. 
The main quantity that discriminates between these two 
different collapse regimes is the flux in the Lyman-Werner 

(LW) bands (11.2-13.6eV). This LW radiation is emitted 

by the first generation of stars and it is convenient to express 
the flux in units of J 21 = 10~®^ ergs“^ cm“^ Hz“^ sr~^ (in 
the following, we use this convention without explicitly writ¬ 
ing J 21 ). The so-called critical value Jcrit sets the threshold 
above which a halo with Tvir > 10^ K can directly collapse 
to a SMBH seed. Below this value, the gas is susceptible 
to fragmentation due to efficient H 2 cooling and the mass 
infall rates towards the centre are generally lower. However, 
the values for Jcrit quoted in the liter ature span several or - 
ders of magnitude from Jcrit = 0.5 (lAear wal et ahl 2015ll 
to as high as Jcrit — 10® llOmukaill200ll : Latif et al.ll2015ll . 
There are several reasons for this large scatter. First of all, 
the value of Jcrit is highly sensitive to the spectral shape of 
the incident radiation field, with softer radiation fields lead- 


ing to significant smaller values of Jcrit (lS hanget_al, 2010|: 

Su ^mura et al.ir2014l : lAgarwal &: Khochfarj|2015l: Latif et ahl 


2 OI 5 II . Secondly, one-zone calculations fe.s. lOmukail l200lf ) 
tend to yield lower values of Jcrit than determinations made 
using 3D numerical simulations. This is a consequence of 
the fact that Jcrit depends to some extent on the details of 
the dynamical evolution of the gas, which are only approxi¬ 
mately captured by one-zone calculations. This dependence 
on the gas dynamics als o leads to Jcrit varying by a factor o f 
a few from halo to halo dShang et aI.ll201^ : lLatif et al ][2Mi). 
Although there seem s not to be one universal value of Jcrit 
(lAgarwal et al.l [20151 ). it is convenient to use this artificial 
threshold as a quantification of the direct collapse scenario 
to test the relevance of different physical processes. Once a 
process significantly affects the value of Jcrit, it is very likely 
that it plays an important role in the formation of SMSs and 
SMBH seeds. 

One of these important processes is H 2 self-shielding 
against LW radiation, which is generally expressed as a sup¬ 
pression factor /sh to the H 2 photodissociation rate (see sec¬ 
tion [2]3]). There are several analytic expressions to calcu¬ 
late /sh as a function of the H 2 column density and the gas 
temperature IPr aine &: Bertoldilll996l : IWolcott-Green et al.1 
l201ll : [Richings et aI.ll2014lL Neglecting self -shielding leads to 
a large change in Jcrit i Shang et al.ll2010r ) and even among 
these analytic fun ctions, the value o f Jcrit varies by an or - 
der of magnitude dLatif et al.1 l2014l : ISugimura et al.l l2014l l , 
which hence shows the importance of a correct treatment of 
this effect. Another challenge is the proper determination of 
the effective H 2 column density for self-shielding, since it is 
either computationally very expensive or not very accurate. 
In this study, we want to test the effect of a more accu¬ 
rate H 2 self-shielding implementation on the direct collapse 
scenario. In contrast to previous studies, we determine the 
column densities self-consistently during the simulation and 
properly account for the Doppler shifts of spectral lines by 
velocity gradients, which reduce the effective column den- 
sity. To do so, we use the treecol algorithm dev e loped by 
IClark et al.1 d2012l ') and extended bv lHa-rtwig et all d2015li to 
calculate a spherical map of column densities around each 
cell. This method is based on the hierarchical tree structure 
used to calculate the gravitational forces between fluid el¬ 
ements in the computational domain and therefore comes 
with only little additional cost. A more detailed description 
of this method is provided in section m 

The paper is organized as follows. In section [2l we 
describe the methods, including initial conditions, chem¬ 
istry network, and the new implementation to determine 
self-shielding. In section [3] we present the results together 
with an analysis of the differences between the self-shielding 
methods and the mass infall rates. We discuss the caveats 
in section [ 4 ] and conclude with a summary in section [5] 


2 METHODOLOGY 

In this section, we present our computational methods. 
First, we explain our initial conditions and refinement strat¬ 
egy. Then, we present the chemical network and our new 
approach for the determination of effective H 2 column den¬ 
sities for self-shielding. 
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Improving H 2 self-shielding 3 


2.1 Initial conditions 

We are interested in the collapse of the subset of metal- 
free haloes that is able to cool by atomic hydrogen. The 
cooling rate of H rises steeply around T ~ 10'* K and ex¬ 
pressed by the virial mass, we are interested in haloes with 
Mvir — IO^Mq. In this study, we focus on the effect of dif¬ 
ferent H 2 self-shielding implementations. Consequently, we 
first run a cosmological dark matter only simulation and se¬ 
lect the first haloes with a mass of ~ 10^ Mq. Under the 
assumption that the value of Jcrit is mainly affected by the 
gas dynamics within the virial radius, this is a representative 
candidate for a metal-free halo that directly collapses to a 
SMBH seed. We assume a flat ACDM Universe and use cos¬ 
molog ical parameters presented by the lPlanck Collaboration! 
(120141 with additional constraints from WMAP polarisa¬ 
tion at low multipoles, high-resolution cosmic microwave 
background data sets, and baryonic acoustic oscillations: 
Ho = 67.8 kms“^ Mpc“^ = h 100 kms”*^ Mpc“^, Ha = 0.69, 

= 0.31, = 0.048, Us = 0.96, erg = 0.83, Une = 0.25. 

We create the initial density field at redshift « = 99 in a peri¬ 
odic box of 1 Mpch“^ comoving with MUSIC dUahn fc Abell 
[20nh . which generates the displacements and velocities fol¬ 
lowing second-order Lagrangian perturbation theory. 

The simulations are pe rformed with t he hydrodynamic 
moving-mesh code AREPO (|SDringe]ll201fll l. which combines 
the advantages of smoothed particle hydrodynamics (SPH) 
techniques and adaptive mesh refinement (AMR) codes. We 
first run a cosmological simulation with ~ 2 x 10^ Voronoi 
cells, which corresponds to a particle mass of rriDM = 

5.1 X IO^Mq. We trace the target halo and a region of 
twice its virial radius to the initial conditions. In a sec¬ 
ond run, we refine this region of interest and also include 
gas, which leads to masses in the highest refined region 
of muM = 100 M© and mgas = 18 Mq . This resolution is 
set to properly resolve the collapse of the gas up to den¬ 
sities of n ~ 10®cm“®. This value is chosen to cover the 
local thermodynamical equilibrium (LTE) of H 2 at around 
n ~ 10^cm“®, because above this value it is much easier 
to collisionally dissociate H 2 than at lower densities. Hence, 
once the gas reaches this value without building up a signif¬ 
icant fraction of H 2 , it is not going to ma nage to do so at 
higher densities either. iRegan et al.l (l2015l l study the effect 
of the dark matter mass resolution that is needed to prop¬ 
erly resolve the collapse of haloes at high redshift. They find 
that for typical collapse scenarios with a moderate LW back¬ 
ground, rriDM = IOOM 0 is a sufficient resolution, whereas 
this minimum mass resolution even decreases for the higher 
LW backgrounds (Jlw = 500) that we want to study here. 
Consequently, our dark matter mass resolution is sufficient 
to properly resolve the collapse. 

The finite box size of our simulations might distort the 
nonlinear effective coupling on the boxscale by not covering 
all relevan t Fourier modes of the power spectrum (see e.g., 
ISetclll999li . However, the choice of our box size is well moti¬ 
vated for a cosmological representative selection of a lO^M© 
halo, because we do not want to draw high-precision cosmo¬ 
logical probes from these simulations, but rather analyse the 
collapse behaviour of one specific halo. Another effect of the 
limited box size is the distortion of the large-scale tidal fields, 
which might affect the angular moment um of the haloes, 
according to tidal torque theory (see e.g.. lFaIl fc Efstathioul 


Il980h . However, the angular momentum budget of haloes 
is dominated by local effects like mergers or the accretion 
of cold gas streams an d only a minor contr i bution comes 
from cosmic tidal fields dPanovich et al.ll2012l : iDubois et al.l 
I 2 OI 4 I : [Uaigle et anbOlSl L In any case, the effect of different 
implementation of H 2 self-shielding on the collapse dynam¬ 
ics should not be affected by the limited box size. 


2.2 Chemistry 


In the following section, we describe our chemical network 
and highlight the most important reactions and rate coef¬ 
ficients. A more extensive discussion of the relevant chem¬ 
ical processes for modelling dire ct collapse wit h a strong 
LW background can be found in iGlove 3 (l2015all^ . We ap¬ 
ply a primordial chemistry network t hat is originally based 
on th e wor k by Glover fc Japps^ (l2007l L I Glover fc Abell 
(l2008tl . and IClark et al.h 2011 ah . Since the deuterium chem¬ 
istry does not affect the direct collapse scenario dGlovetl 
l2015ali . we only follow explicitly the evolution of H, He, H 2 , 
H"*", H“, , He^, He^’*', and e~. iGloveil d2015ai i identi¬ 

fied a minimal subset of reactions that must be included in 
the chemical model in order to determine Jcrit accurately. 
We have made sure to include all of these reactions in our 
chemical network. Full details regardin g our choice of reac ¬ 
tion rate coefficien ts can be found in iGlark et al.l d2011al i 
andiCT ove 3 d2015al l^. 

The collapse dynamics depends strongly on the abun¬ 
dance of molecular hydrogen, which is the dominant coolant 
for temperatures below ~ 10"^ K. Molecular hydrogen is 
mainly formed via the two-step process: 


H -f e —>• H -I- 7 
H-hH" H 2 -f e". 


( 1 ) 

( 2 ) 


and is primarily destroyed either by collisions with hydrogen 
atoms 


H 2 -f H ^ H -f H -f H, (3) 

or by the so-called Sol omon process (iField et al.l Il966l : 
IStecher fc Williamslll967h : 

H 2 - 1 - 7 ^ H -h H, (4) 

where a LW photon photodissociates the molecule by excit¬ 
ing it from the electric ground state into an excited electronic 
state. In ~ 15% of the cases, the electrons do not decay 
into a bound state, but into the vibrational continuum of 
the ground state and thereby dissociate the molecule. For 
the Solomon process, we assume a spectrum of a blackbody 
with an equivalent temperature of T^ad = 10® K (T5) that 
is cut oft above 13.6 eV, because photons with higher ener¬ 
gies are absorbed by the intergalactic medium. This choice 
represents the case in which the spectra are dominated by 
Pop HI stars. The second generation of stars is believed to 
be less massive and can be approxi mated by a blackbody 
spectrum with 10"^ K < Tra d < 10® K (ISugimura et al.ir2014l : 
lAgarwal fc Khochfarj|201^ . It is important to note that the 
value of Jcrit depends on the choice of the spectrum, be¬ 
cause whereas a T5 spectrum mainly photodissociates the 
H 2 directly (equation , a spectrum with a cooler effective 
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temperature dominantly prohibits H 2 formation by photode¬ 
tachment of H“ via 


H -f 7 —^ H -I- e 


( 5 ) 


jLatif_et_aJj _ 1201^: ISueimura et ap I20l4 

lAgarwal fc Khoch^ 20151 : lAgarwal et al. 201^. Since 
we want to focus on the effects of H 2 self-shielding, we 
will only consider the T5 spectrum, where this effect plays 
a dominant role. We also include d i ssocia tive tunnelling, 
which is discussed in iMartin et al.l (Il996fl . This process 
significantly contributes to the total collisional dissociation 
rate of H 2 and is therefore necessary for a proper treatme nt 


_ L p n __ 

of primordial gas physics (iLatif et al.ll2014l : lGloverll2015al l. 


2.3 H 2 self-shielding 

Since one photon of the external radiation field can only 
photodissociate one H 2 molecule, a large column density of 
molecular hydrogen can protect the inner regions against 
photodissociation. This process is known as self-shielding. 
For the photodissociatio n of H 2 , we use the rate coefficient 
dGlover fc Japps^l2007l ~l 

= 1.38X ( 6 ) 

J 21 

where the factor fsh ^ 1 accounts for the effect of H 2 self¬ 
shielding with fsh = 1 in the optically thin limit. This rate 
coefficient corresponds to a normalisation of the radiation 
field of Jlw = 1 at 12.87 eV, the middle of the LW bands. 
The exact treatment of H 2 self-shielding requires a full ra¬ 
diative scheme with line transfer of all the important lines in 
the LW bands and is therefore prohibitively expensive. How¬ 
ever, the shielding factor can be approximately expressed as 
a function of the H 2 column density and the gas temper¬ 
ature. The latter enters because of the thermal broaden¬ 
ing of spectral lines and due to the temperature-dependent 
excitation of different ro tational levels of the H 2 molecule. 
iDraine fc Bertoldil (Il996l l study the structure of stationary 
photodissociation fronts and propose a self shielding factor 
of the form 


the function by IPraine fc Bertol dil d 199 611. A n other func¬ 
tional form was proposed by Richings et al.l (|2014| ~) who 
model shielding against UV radiation in the diffuse inter¬ 
stellar medium. They also include the effect of turbulent 
gas velocities, but they derive the self-shielding factor for a 
one-dimensional plane-parallel slab of gas. A lthough the fit¬ 
ting function by I Wolcott-Green et al.l (1201 ll i was derived in 
a three-dimensional simulation, it is also based on a static 
slab of gas. 

In a more realistic scenario, however, we are interested 
in a collapsing cloud, where the relative velocities between 
infalling gas particles Doppler-shift the spectral lines. Due 
to this effect, an H 2 molecule can only shield other H 2 
molecules whose relative velocity is smaller than its thermal 
velocity. Otherwise, the spectral lines are shifted too far and 
H 2 molecules do not contribute to the effective colum n den¬ 
sity. To account for this effect, iHartwig et al.l ll2015ll have 
implemented a new method for the determination of effec¬ 
tive column densities in three-dimensional simul ations. This 
metho d is based on the treecol algorithm bv IClark et al.l 
l| 2012 h . which directly sums up the individual mass contribu¬ 
tions for the column density of each fluid element. With this 
information, we create spherical maps of the column density 
around each Voronoi cel l, with 48 equal-are a pixels based on 
the HEALPIX algorithm dGorski et al.l[2005l l. The number of 
48 pixels is motivated by healpix, which divides the sphere 
into 12 equal-area pixels, which can be subdivided int o 2^ 
pixels each. Based on the work by iGlark et al.l (l2012fl . we 
chose N = 2, since this value provides a sufficient angular 
resolution for most astrophysical applications. 

Two characteristics of the code make treecol highly 
useful for our purpose. First, it uses the tree structure, which 
is already present in the code to determine the gravitational 
force. Hence, the determination of column densities comes 
with only little additional computational cost. Secondly, we 
can directly compare the relative velocities of the particles Vr 
that possibly contribute to the column density to the ther¬ 
mal gas velocity uth of the particle for which we want to ca l- 
culate the column density. Following IHartwig et al.l ll2015ll . 
a particle only contributes to the effective column density, if 


^ _ 0.965 , 0.035 

“ [l + x/bhY ^ (l-f 2;)0'5 

where x = AH 2/5 x 10^'^cm~^ is the H 2 column density, 
65 = ti/10®cms“^ is the scaled Doppler parameter of the 
molecular hydrogen, and a = 2. This functional form ac¬ 
counts for the effect of line overlap and has been applied in 


manv studies (iGlover & JaDDsenI l2007l: 

Whalen & Norman 

2008 

: Gnedin et alJ l2009l: Shang et al.l 

2 OIOI: Glover et al. 

2010 

: Christensen et al.l 20121: Krumholz 

2OI2I1. 


Wolcott-Green et al.l (201lf) model 

the photodissocia- 


tion and self-shielding of H 2 in protogalaxies with three- 
dimensional simulations, based on post-p rocessing their out¬ 
put da ta. They found that the formula bv iDraine fc Bertolda 
(Il996ll is only valid for cold or low density gas, in which only 
the lowest rotational states of H 2 are populated and propose 
the modification a = 1.1 to equation 0 , which provides a 
better fit in dense gas for all relevant temperatures. A com¬ 
paris on of both functions dLatif et al.ll20l4 ISugimura et al.l 
2014|) shows that th e appli cation of the newer formula by 
Wolcott-Green et al.l (I 2 OIIII yields values of Jcrit that are 
up to an order of magnitude lower than those derived with 


Vr < 1.694uth, ( 8 ) 

where wth is the thermal gas velocity and the numerical fac- 
tor is an exten sion of the so-called Sobolev approximation 
llSobolevlll960h and takes the true line profile into account. 
Based on this criterion, we determine the effective column 
densities and calculate the shielding factors separately for 
all pixels with equat i on (TTll and the exponent a = 1.1 by 
IWolcott-Green et al.l (1201^ . The final shielding factor is the 
mean of 48 directional-dependent factors. 

With this new approach, w e can use formula 0 with 
the IWolcott-Green et al.] 1I2OIIII exponent a — 1 . 1 , which 
was derived for static gas and extend it to collapsing gas 
clouds by defining the effective H 2 column densities based 
on the relative gas velocities. Since our approach automati¬ 
cally accounts for turbulent motions on scales above the spa¬ 
tial resolution, the Doppler parameter 65 in formula 0 does 
only include the thermal broadening. This approach cannot 
only be used for H 2 self-shielding, but also for many other 
radiative transfer processes that rely on the determination 
of column densities. The only requirement is that the reso¬ 
lution elements are stored in a tree-like structure, which is 


(i+*r 


1180 


, ( 7 ) 
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already the case in most codes that include self-gravity. This 
method for the determination of e ffective column densi ties 
is tested and explained in detail in iHartwig et ahl ll2015ll . 

Our method of computing effective H 2 column densi¬ 
ties is valid as long as the main contribution comes from 
the core of individual lines that shield themselves. At high 
column densities, however, the Lorentzian contribution to 
the line profile becomes important and the corresponding 
damping wings should be taken into account for the de¬ 
termination of self-shielding . While this effect is neg ligible 
at small column densities, iGnedin fc Drain'3 (12014 1 show 
that it should be taken into account for H 2 column densities 
of Nb _2 > 10^^ cm“^. At these high column densities, the 
Doppler-shifts induced by relative velocities become less im¬ 
portant and eventually the total column density contributes 
to self-shielding. A more detailed analysis of this effect is 
given in section [SA] where we show that a correct treatment 
of the overlap of these damping wings changes the value of 
the self-shielding factor by less than 5% and has no influence 
on the determination of Jcrit. In addition, the use of effective 
column densities is computationally more efficient, since the 
velocity criterion imposed by equation (|8} limits the amount 
of fluid elements that have to be projected. Therefore, we 
employ this method in our numerical simulation s. _ 

We also note that I Safranek-Shrader et al.l (l2012li use 
a non-local approach for the determination of H 2 column 
densities. They study the influence of LW radiation on star 
formation in the first galaxies and approximate the column 
density in the following way. For each computational cell, 
they calculate the column densities in six directions parallel 
to the coordinate axis. The smallest of these column densi¬ 
ties is then used to calculate the self-shielding factor. They 
show that this is already an improvement, compared to lo¬ 
cal estimations of the column density. Similar techniques 
have also been used to study the ef fects of H 2 self-shielding 
in giant molecular clouds Isee e.g. iNelson fc Langen 119971 : 
iGlover fc Mac L^l20078] fal. 

Most other previous simulations use an approximation 
for the H 2 column density, based on a characteristic length 
scale Lchar, and the assumption that the H 2 density is con¬ 
stant within this length and negligible beyond it. The col¬ 
umn density is then given by Ahs = UHjAchar, where uhs 
is the local number density of molecular hydrogen. Assum- 


ing that the effe 
simulations (e.j 

5ct of self-shielding occurs only locally, many 

r. IShanff et al.l 12010 

: IVan Borm et al.l l2014l: 

Sueimura et al. 

2014 

: iLatif et al.l 

2OI5I: iLatif & Volonteri 

2 OI 5 I: lGloverll2015allb 

: lAffarwal et al.||2015|'l use = Lt 


with the Jeans length 


Aj 


IS/cbT 
4'KGfj,mpp ’ 


(9) 


where T is the temperature, p the mean molecular weight 
and p the gas density. Since this is a widely used approxima¬ 
tion, we will compare our results obtained with TREECOLto 
this formula. 


3 RESULTS 

In order to increase the statistical significance of our results, 
we create four independent sets of cosmological initial con¬ 
ditions. The side length of each of these boxes is lMpc/i“^, 



-1000 -500 0 500 1000 


pc 



Figure 1. Maps of the average number density of hydrogen nuclei 
along the line of sight for halo C with the TREECOL approach at 
the moment of collapse {z ~ 15.1) at a scale of 2000 pc (top) and 
20pc (bottom). The background flux is Jew = 10^ < Jcriti and 
we can clearly see the formation of clumps in the central region. 


and we select one halo per box. As described in section E7I1 
the region of interest is refined for each box and we resimu¬ 
late this set of cosmological zoom-in simulations. At redshift 
a = 30, we switch on the photodissociating background. In 
reality, the LW radiation increases with cosmic time and the 
time of its onset also influen ces the collapse of primordial 
gas clouds (IVisbal et al.ll2of3 ). We use this simplification of 
an instantaneous onset to be able to focus on the different 
implementations of the H 2 self-shielding and make it com¬ 
parable to most previous works in this field. As long as the 
H 2 abundance has enough time to reach an equilibrium, be¬ 
fore the halo of interest starts with the run-away collapse, 
the results are unaffected by the choice of this redshift. This 
criterion is fulfilled in all our simulations. 

Altogether, we study the collapse of four different haloes 
(A, B, G, and D) with two different methods (treecol and 
Jeans approximation) for determining the column density 
and several different strengths of the LW background per 
halo. All these runs are completely independent, and the 
column densities are calculated self-consistently during the 
simulation. As an example, the structure of halo G and its 
central region is illustrated in Fig. [T] We can see that the 
halo is embedded in the cosmic web and is fed by several 
gas streams. The central region shows a lot of substructure 
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and gas clumps, which indicate ongoing fragmentation. In 
this case, the LW radiation is not strong enough to prevent 
efficient H 2 cooling and the gas can locally contract before 
the cloud globally collapses. 


3.1 Determination of Jcut 

The value of Jcrit sets the threshold between the two differ¬ 
ent collapse regimes. Above this value, the H 2 fraction re¬ 
mains low, the temperature stays around 10'* K during the 
collapse, and only one central density peak forms. Below 
this value, the photo-dissociating background is not strong 
enough and H 2 line emission cools the collapsing gas to a few 


hundred K. This collapse tyi 

sically results in several frag- 

ments (IClark et al. 

2 OIIJ 1 J: 

Greif et al.ll201ll: IStacv et al.l 

I 2 OI 2 I: iHirano et al.l 

2014|). In order to discriminate between 


these two scenarios, we have to analyse the temperature evo¬ 
lution during the collapse. For two typical cases, the phase 
diagrams are given in Fig. (2] During the virialisation of the 
halo, the gas shock heats to the virial temperature of around 
10^ K. The gas contracts further and remains at this temper¬ 
ature due to cooling by atomic hydrogen. With increasing 
density, also the column density of H 2 increases and the self¬ 
shielding against the external radiation becomes more effi¬ 
cient. As discussed above, H 2 reaches LTE at n ~ 10"* cm“®, 
and if the collapse is still isothermal up to these densi¬ 
ties, it will proceed isothermally. If the LW background is 
not strong enough, the H 2 fraction can increase, which in 
turn increases the column densities and hence leads to a 
more efficient self-shielding, which consequently increases 
the H 2 fraction even further. This runaway production of 
H 2 enables a clear distinction between the two different col¬ 
lapse regimes. An H 2 fraction of ~ 10“^ is sufficient to cool 
the gas to temperatures of a few hundred K and hence to 
induce gas fragmentation. The individual clumps that form 
in this latter scenario can also be seen in the phase diagram 
as stripes in the cold, high-density regime, indicating that 
their thermal evolution is decoupled from one another. 

In order to find the values of Jcrit, we have to com¬ 
pare the phase diagrams for different values of the LW back¬ 
ground. To be able to do so, we bin and plot the data in log¬ 
arithmic density space. An iterative algorithm merges and 
splits these density bins until each bin contains roughly the 
same number of Voronoi cells. This procedure ensures a sta¬ 
tistically representative binning, but can have the side ef¬ 
fect that the actual peak density in the simulation is not 
the same as the one in the binned data. The correspond¬ 
ing plots for the four haloes and the two different methods 
are shown in Fig. [21 We test and display the background 
strengths of Jlw = 10 ®, 3 x 10®, 10"* for all haloes and then 
successively bracket the actual value of Jcrit. We stop the 
simulations after the first snapshot with a peak density of 
n ^ 10® cm~® and compare the temperatures in the density 
regime above n ^ 10^ cm“®. If the temperature falls below 
6000 K, the collapse is regarded as non-isothermal. With this 
method, we find the lowest value Ji for which the collapse 
is still isothermal and the highest value J2 for which H2 can 
efficiently cool the gas. Due to a limited number of possible 
realisations, we define the final critical value as the geomet¬ 
rical mean between these two values 

Jcrit = VJi J 2 . (10) 
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Figure 2. Temperature as a function of number density of hy¬ 
drogen nuclei colour-coded by the H 2 fraction for halo C using 
TREECOLwith Jew = 10® (top) and Jew = 10"* (bottom). To 
better illustrate the relevant difference at higher H 2 fractions, we 
artificially set a lower limit of /hj = 10“t0 in these plots. We 
can clearly see the different collapse behaviours depending on the 
strength of the LW background. With a high JeWi the gas re¬ 
mains hot around 10^ K with /jjj S) 10“^. For a lower strength 
of the photodissociating background, the fraction of molecular 
hydrogen rises up to /jjj — 10“® and the gas cools down to sev¬ 
eral hundred K. In the latter case, small clumps decouple from 
the global thermal evolution and we can see their footprints as 
stripe-like structures in the cold high-density gas. 


Due to this finite number of tested Jew, the final val¬ 
ues have an uncertainty of ~ 10%. The resulting values, 
the virial masses of the haloes, and the collapse redshifts 
are compared in Table [T] First of all, we directly see that 
Jcrit is about a factor of two lower in the runs where we 
use the treecol method to calculate the column densities. 
The reasons for this effect will be discussed in detail be¬ 
low. However, already the results based on the commonly 
used Jeans approximation are lower than found in previ¬ 
ous studieS;_Thereare t wo main reason s for this. First, 
e.g. lLatif et al.l (l2014l ~l and lSugi mura et al.l (l20 14h sho w that 
the self-shielding function bv IWolcott-Green et al.l (1201 il l 
yields values of Jcrit that are up to an order of magni- 
tu de lower, compared to th ose derived with the function 
by iDraine fc Bertoidil lll996fl . Secondly, the ENZO chemical 
model, which was used by many studies in this field (e.g. 
IShang et al.l l20ld : I Wolcott-Green et al.l I 2 OIII : [Regan et al.l 
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Figure 3. Temperature as a function of the number density of 
hydrogen nuclei for the four haloes and the two different meth¬ 
ods for determining the H 2 column density. The curves represent 
several realisations with different LW backgrounds, where the 
long-dashed lines represent the isothermal collapse and the short- 
dashed lines the collapse with efficient H 2 cooling. From these 
plots, we can read off the critical value for the isothermal col¬ 
lapse, which is systematically lower for the runs based on the 
TREECOL method. 


l2014l : [Latif et alj 2014'), tends to overestimate Jait by about 
a factor of two jGloveill2015al l. Hence, our results based on 
the Jeans approximation are in compliance with previous 
studies. 


3.2 Differences in the H 2 self-shielding 

In order to understand the differences induced by the new 
TREECOL approach, we have to compare the column den¬ 
sities and corresponding self-shielding factors for the two 
methods. First of all, we summarise the most important fea- 
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halo 

Jcrit (Jeans) 

Jcrit (tREECOL) 

Mvir/ 1 O®'M 0 

■^coll 

A 

3500 

1700 

1.8 

17.9 

B 

5500 

2200 

1.2 

14.4 

C 

3500 

2200 

1.7 

15.1 

D 

3500 

1700 

1.1 

12.9 


Table 1. Critical values Jcrit of the LW background for the four 
different haloes and the two different column density approaches. 
The Jcrit determined with the TREECOL method is smaller by 
about a factor of two in all haloes. The halo-to-halo variance 
of this value is small and the collapse redshifts are distributed in 
a reasonable range. We also list the virial mass and the collapse 
redshift, which indicates the time when the halos first reach a 
density of n ^ 10® cm“®. 


tures of this new method to understand these results. With 
TREECOL , we create a spherical grid with 48 pixels around 
each Voronoi cell and project the column densities on to this 
grid. However, we do not calculate the total column densi¬ 
ties, because gas can only contribute to the self-shielding if 
its relative velocity is smaller than ~ 1.7 times the thermal 
velocity (equation |8]). Based on this criterion we determine 
the spherical maps of effective column densities, which rep¬ 
resents the spatial distribution of the self-shielding gas. The 
self-shielding factors are then calculated based on the col¬ 
umn density for each of the pixels separately and then av¬ 
eraged over the 48 directions. This procedure is physically 
motivated, because the product Jlw/sL in equation m rep¬ 
resents an effective photodissociating flux and we simply 
average these fluxes over 48 different directions. 

The importance of this point becomes clear, if we anal¬ 
yse the directional dependence of the column densities in 
Fig. |4] We see that the column density distribution around 
one Voronoi cell is generally highly asymmetric and domi¬ 
nated by the contributions from one direction. In the case of 
a collapsing halo it is certainly the central high-density peak 
that yields the strongest contribution to the column density. 
From this direction, obviously, we should not expect any 
photodissociating radiation. A different way of presenting 
this important directional dependence of the self-shielding 
factor is given in Fig. [5] Averaging the column densities first 
and then calculating the shielding factor based on the one 
mean column density is highly biased by the contribution 
of one dominating direction. Consequently, it is important 
to properly average the effective photodissociating fluxes, 
as we do in our simulation. In contrast, a local approxima¬ 
tion that only yields one column density might be affected 
by the central density peak and consequently underestimate 
the self-shielding factor. 

Now, we want to compare directly the results for the 
column density and the shielding factor based on the two 
implementations. In Fig. we display the column densities 
of the simulations with the Jeans approximation and the 
TREECOL approach. These column densities are calculated 
self-consistently during the run and are therefore not di¬ 
rectly comparable. However, we already see that the column 
densities in the isothermal case remain under 10®® cm“® and 
that, for smaller values of Jlw, the column density increases 
already for lower densities. To be able to directly compare 
the effect of the two different approaches, we use one snap- 
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Figure 4. Plots of the maximal column density for each Voronoi 
cell A^h 2 ( 1 ) a function of the mean H 2 column density aver¬ 
aged over the remaining 47 pixels A^h2((2 : 48)), colour-coded 
by the radius. In the top panel, we show the results for halo C 
and a flux of Jlw = 10 ^ < Jcrit ^-nd in the lower panel for the 
same halo but with Jlw = 10"^ > Jcrit- This plot illustrates the 
huge directional dependence of the column densities, because for 
one Voronoi cell, the column densities in different directions vary 
by up to ten orders of magnitude. Interestingly, we can also see 
the more spherically symmetric collapse structure in the isother¬ 
mal case (lower panel), because the plot converges towards the 
diagonal for small radii, indicating that the column densities are 
distributed isotropically close to the centre. Whereas in the up¬ 
per panel, the angular distribution of column densities remains 
anisotropic even in the central region. 


shot of the simulations based on treecol and determine the 
corresponding Jeans column densities by post-processing the 
data. The comparison of these column densities is given in 
Fig. [7| In the isothermal case, the column densities remain 
below 10^® cm“^ and in this regime, the Jeans approxima¬ 
tion overestimates the column densities by approximately 
one order of magnitude. Consequently, the shielding in this 
regime is more efficient with the Jeans approximation, and 
we need a higher LW background to obtain an isothermal 
collapse. For column densities above 10^®cm“^, the Jeans 
approximation underestimates the values. We find that these 
regions, where the Jeans approximation underpredicts the 
column densities are the lower density regions between the 
clumps (compare e.g. with the lower panel in Fig. In 
this low-density environment, the Jeans approximation sees 


Figure 5. Self-shielding factor for Halo C with the 
TREECOL approximation and Jlw = 10^ < Jcrit? colour-coded by 
the gas temperature. This plot illustrates the directional depen¬ 
dence of the self-shielding factor. In our simulation, we determine 
individually the self-shielding factors for the 48 different direc¬ 
tions and average them afterwards (horizontal axis of this plot). 
On the vertical axis we see the self-shielding factor for exactly 
the same simulation output, but here, we first average the 48 dif¬ 
ferent column densities for each fluid element and calculate the 
self-shielding factor based on this one averaged column density. 
The proper direction-dependent treatment of the H 2 self-shielding 
generally yields a less efficient shielding against LW radiation. 
This discrepancy is smaller for temperatures of a few hundred K. 


Jeans TreeCol 



log(nH / cm"^) 109(0^/001"^) 

Figure 6. Molecular hydrogen column density as a function of 
the density of hydrogen nuclei for halo C. For the TREECOL runs, 
we plot the median column density over the 48 pixels in order 
not to be biased by one dominating direction. The values of the 
column density are calculated self-consistently during the simula¬ 
tion. Hence, they are not directly comparable for the two different 
approaches, because the structure of the cloud might be differ¬ 
ent. However, we can already see important similarities such as a 
threshold column density of about 10^® cm~^ above which strong 
self-shielding enables the formation of sufficient H 2 to cool the 
cloud efficiently. 


only the local gas conditions and predicts a rather small col¬ 
umn density, whereas treecol is able to capture the nearby 
high-density clumps with high H 2 fractions. Consequently, 
TREECOL (correctly) yields higher column densities in this 
regime. However, this underestimation does not affect the 
value of Jcrit, because we only get H2 column densities sig¬ 
nificantly above 10 ^® cm“^ in haloes where Jlw < -Jcrit and 
where the gas can undergo runaway cooling. Phrased dif¬ 
ferently, this underestimation by the Jeans approximation 
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Figure 7. Direct comparison of the H 2 column densities for the 
two different approaches. The data presented here is based on an 
output of halo C simulated with TREECOL and the corresponding 
column densities based on the Jeans approximation are calculated 
based on the output file. Again, we use the median over the 48 
pixels for the column densities based on TREECOL . For column 
densities below 10^®cm“^, the Jeans approximation yields col¬ 
umn densities that are higher by about one order of magnitude. 


Figure 8. Direct comparison of the self-shielding factors for 
the two different approaches. The data presented here is 
based on haloes simulated with TREECOL and the correspond¬ 
ing Jeans column densities are calculated based on the out¬ 
put files. In the isothermal regime (long-dashed lines), the 
TREECOL approximation yields higher self-shielding factors and 
therefore provides a less efficient shielding against the external 
photodissociating radiation field. 


at higher column densities is a consequence of the runaway 
H 2 cooling and the subsequent fragmentation and not its 
trigger. 

There are two obvious reasons why TREECOL yields 
lower effective column densities. First of all, it takes into ac¬ 
count the three-dimensional distribution of the matter and 
the declining density radially outwards, whereas the Jeans 
approximation assumes a constant density within one Jeans 
length and hence generally overestimates the gas number 
density. Secondly, only fluid elements within a certain ve¬ 
locity range contribute to the effective column density in 
TRE ECOL , which again reduces t he value of the column den¬ 
sity. IWolcott-Green et al.l ll201ll ) have already pointed out 
that the Jeans approximation generally overestimates the 
column densities in a static, isothermal slab of gas. Our 
treatment of the relative velocities further increases this ef¬ 
fect. 

We can see the same trend for the self-shielding fac¬ 
tors in Fig. [S] In the case of an isothermal collapse, the 
TREECOL methods yields always higher values for the shield¬ 
ing factor and therefore a less efficient shielding against the 
LW background. Consequently, a smaller value of Jcrit pro¬ 
vides already enough photodissociating radiation to keep the 
collapse isothermal. 

3.3 Impossibility of a simple correction factor 

Although TREECOL is computationally efficient, it would be 
good if there is a simple correction formula to the Jeans 
approximation that is able to reproduce the self-shielding 


results obtained from TREECOL at least to some degree. The 
aim would be to find an updated version of equation 0 that 
uses the column densities A^H 2 = based on equation 

®. As a start, we can impose the velocity criterion from 
equation ([ 8 ]) on the matter included in the calculation of the 
effective H 2 column density. This changes the dependence 
of the self-shielding factor on the thermal velocity, which 
is expressed by the exponent a in equation ©. Hence, we 
try to fit the results obtained with TREECOL with the free 
parameter a* but with the column densities determined with 
the Jeans approximation as an input: 

/sh (A^H2 jTreeCoOfi — 1.1 — ,Jeans)^* (11) 

However, this exponent varies with a huge scatter between 
the different haloes and with the various strengths of the 
background radiation. For example, for halo C we find a* = 
0.6 for Jew = 2400 and a* = 1.2 for Jew = 4000, although 
both collapse isothermally. Most other values are distributed 
in the range 0.7 < a* < 1.1, which illustrates that it is not 
possible to find a simple correction factor to reproduce the 
results by TREECOL . This is mainly based on the fact that a 
proper determination of the self-shielding factor has to take 
into account the three-dimensional structure of the density 
and velocity field, as we have seen before. Since the Jeans 
approximation is only based on local quantities, it cannot 
capture this structure and consequently fails at reproducing 
the self-shielding factors. 


© 2015 RAS, MNRAS OOO.ITlfTil 
























10 


T. Hartwig et al. 



4.0 


3.5 




3.0 05 


2.5 

2.0 


Figure 9. Total H 2 column density as a function of the effective 
H 2 column density colour-coded by the gas temperature for halo 
C with Jlw = 1000 < Jcrit- For’ tho effective column density, 
we impose the velocity criterion of equation [8| and only matter 
fulfilling this criterion is included. The excess at A^H 2 ,treecol — 
10^®cm~^ is created by the hot gas, which falls on to the cold 
and dense gas clumps and has therefore a high relative velocity. 
Hence, the velocity criterion excludes a lot of matter, and the 
total column density is significantly higher in this regime. 


3.4 Effect of damping wings 


So far, we assumed that only matter in a certain veloc¬ 
ity range contributes to the self-shielding of H 2 , because 
for larger relative velocities the spectral lines are Doppler- 
shifted too far from the core. As already mentioned in 
section El this picture changes for higher column den¬ 
sities, where the contribution of the damping wings be- 
comes important and different lines can self-shield e ach other 
l|Black fc Dalgarn3ll977l : IPraine fc BertoldilllQQtt . This ef¬ 
fective broadening of spectral lines makes relative velocities 
less important and the total H 2 column density then con- 
tri butes to self-shieldi n g. Th e original self-shielding formula 
bv IPraine fc Bertoldil lIlQQtt . equation E in section El al¬ 
ready included these two contributions, where the first term 
corresponds to shielding from the cores of individual lines, 
and the second term represents the effect of overlapping 
dam ping wings, which are n ot affected by relative veloci¬ 
ties. [CnSl^TiDrai^ ||2014 1 propose a correction to the 
self-shielding factor of the form 


fsh 


0.965 0.035 


(1 + X2)°-" ' 
1180 


, ( 12 ) 


where X\ = AiH2.treecoi/5 X lO’^'^cm”^ is the effective H 2 col¬ 
umn density and X 2 = VHo.totai/5 x 10 ^^cm ~^ is the total 
H 2 column density. iGnedin fc Drain'3 ([20^ ) use the expo¬ 
nent a = 2 based on the original work by Draine fc Bertoldil 
lll996h. whereas we apply t he more recent value a = 1.1 by 
I Wolcott-Green et al.l (1201 il l. 

We investigate the impact of the damping wings by de¬ 
termining the total H 2 column densities. To do so, we post 
process the snapshots of halos G with treecol and compare 
them to the effective column densities (Fig. [9]). The total col¬ 
umn density is higher by up to several orders of magnitude. 
Especially around Vh2, treecol — 10'^®cm~^, we see that the 
total column density is much larger than the effective one. 


This excess is created by the hot gas that falls towards the 
centre and has therefore a high relative velocity with respect 
to these central cold and dense gas clumps. Consequently, 
the velocity criterion excludes a large contribution, and the 
total column density is significantly higher than the effective 
one in this regime. In any case, this comparison strengthens 
the importance and illustrates the influence of the velocity 
criterion on the determination of the column densities for 
self-shielding. 

We now use these two column densities and determine 
the self-shielding factors with equation m and compare 
them to the previously used values based on equation (E. In 
the case of Jlw = 10“^ > Jcrit, the effective and the total col¬ 
umn densities remain below ~ 10^®cm“^, and the contribu¬ 
tion of damping wings is negligibly small: the mean deviation 
between the self-shielding factors is 0.24% with a maximal 
discrepancy of 0.60%. This is due to the fact that the main 
contribution to the self-shielding in this regime comes from 
the cores of individual lines and hence, from the first term in 
equation (inj. For the case with Jlw = 10® < Jcrit, we have 
a higher H 2 abundance and consequently higher H 2 column 
densities. Here, the self-shielding factors differ on average 
by 2.2% with a maximal deviation of 5.3%. Hence, we con¬ 
clude that the contribution of the damping wings to the 
self-shielding is negligibly small in our scenario, compared 
to other effects and approximations. In particular, it seems 
not to affect the determination of Jcrit, because this effect 
only becomes important, once the runaway H 2 production 
has already set in. 


3.5 Mass infall rate 


The main quantity that determines if a SMS and hence 
a seed of a SMBH forms is the mass infall rate Afin , 
which has to be above Min S*. O.lMp) llBegelmanI 20101: 
HosokwaetaL 20131: Schleicher et al.l 120131 : iLatif et ^ 


2013bl : Ferrara et al.ll20i3) . An isothermally collapsing cloud 
at T ~ 10"* K is a sufficient criterion to provide the neces¬ 
sary mass infall rates, but is this criterion also necessary? 
To study this question, and to see if the H 2 self-shielding 
implementation induces any difference, we analyse the mass 
infall rate for the different haloes and methods. Assuming 
a spherically symmetric cloud, we determine the mass infall 
rate by 


Min = 47rr®|ur|p, 


(13) 


where r is the distance from the densest point, Vr is the ra¬ 
dial velocity, and p is the gas density. We average the radial 
velocity and the gas density over spherical shells. However, 
we note that the mass infall rate is generally not spherically 
symmetric and neither is it constant in time. We therefore 
take the mean over the last ~ 10® yr of evolution to ob¬ 
tain a reasonably smooth and well-defined value (Fig. llOl) . 
The mass infall rates in the isothermal collapse scenarios 
are in the range 0.1MQyr“® < Min < 1 MQyr“® and tend 
to be higher than those in the cases with Jlw < Jcrit. Espe¬ 
cially in the central ~ iQ pc, only the isothermal clouds are 
able to provide mass infall rates above 0.1MQyr“®, which 
are necessary for the formation of SMBH seeds. As we have 
seen before, the self-shielding method influences the value of 
Jcrit, but if we now compare the regimes below and above 
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Figure 10. Radial profiles of the mass infall rates averaged over 
the last ~ 10 ® yr of the collapse. Simulations with Jlw < ■.^crit 
are shown with short-dashed lines and those with Jlw > •.^crit 
are shown with long-dashed lines. The black line at Min = 
0.1 MQyr“^ represents the theoretical threshold above which the 
formation of a SMBH seed is possible. Generally, the mass infall 
rates are higher in the case of an isothermal collapse and espe¬ 
cially in the central regions, the mass infall rates fall significantly 
below O.lM 0 yr“^ for Jlw < •.^crit- 


Jcrit separately for both approaches, we do not see any clear 
systematic difference that is induced by the choice of the 
self-shielding method. 

Assuming that approximately one Jeans mass per free 
fall time falls towards the central region, the mass accretion 
rate is given by 

3 

M = ^ oc (14) 

Gr 

where Cs is the sound speed and G the gravitational con¬ 
stant. Consequently, the accretion rate is higher for hotter 



Figure 11. Map of the average number density of hydrogen nuclei 
along the line of sight for halo C with the TREECOL approach at the 
moment of collapse (z ~ 15.1). The background flux is Jlw = 3 x 
10 ® > Jcrit and in contrast to Fig. [J we can see no indication for 
gas fragmentation but a rather smooth, approximately spherically 
symmetric accretion towards the centre. 


gas, such as in the isothermally collapsing cloud that re¬ 
tains temperatures of ~ lO'^ K during the collapse. More¬ 
over, the Jeans mass drops with the temperature and the 
gas is more susceptible to fragmentation in the case of effi¬ 
cient H 2 cooling (see e.g. IClark et al.l[2011al [bh. Once cooling 
becomes efficient and the gas temperature falls, the cooling 
time becomes shorter than the local freefall time and the gas 
can locally contract, before the cloud globally collapses. The 
resulting clumpy structure can be seen in Fig.[Tl whereas the 
spatial structure of the isothermal collapse is illustrated in 
Fig. 111! This smooth, almost spherically symmetric struc¬ 
ture without signs of fragmentation enables higher gas infall 
rates, which leads to the formation of a SMS which then 
collapses to a SMBH seed. However, one should keep in 
mind that the different mass infall rates for scenarios be¬ 
low and above Jcrit are just a trend and also the threshold 
of 0.1MQyr“® is only a rough estimator with other pro¬ 
posed values between 0.01 and 1 M^yr"® llBegelmaniBoid: 
Hbsokawa et al.l I2OI2I : ISchleicher et al.l I2OI3I : Ferrara et al.l 


2 OI 4 I) . We find Min < O.lMQyr ® also for the isothermal 


collapse and vice versa. A value of 0.6 M 0 yr“^ seems to dis¬ 
criminate between the two collapse regimes: all isothermally 
collapsing clouds yield accretion rates above this value, but 
only one out of eight clouds with Jlw < Jcrit reaches this 
accretion rate. 

Recently. iLatif fc Volont^ (l2015l ) study the infall rates 
in atomic cooling haloes in greater detail. As in our study, 
they find that it is not always necessary to completely sup¬ 
press H 2 formation to obtain sufficiently large infall rates to 
form a SMS. Moreover, they detect a rotationally supported 
structure in the central parsec, but this rotational support 
does not halt the collapse and still enables infall rates of 
~ O.lM 0 yr“®. For a more detai led discussion of t h is top ic, 
we refer the interested reader to lLatif fc Volonteril (l2015l ). 
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4 CAVEATS 

The attempt to find one universal value of Jcrit is rather 
artificial, because the relevant physical processes are too 
complex to be summarised in one simple number that de- 
cides whether we fo rm a SMBH seed or not. Recently, 
lAgarwal et ahl (l2015l l study the value of Jcrit using one-zone 
models with a more realistic spectral energy distribution for 
the external radiation and show that Jcrit is not one fixed 
value, but rather a range spanning more than three orders 
of magnitude. Moreover, we can also have sufficiently high 
mass infall rates to form a SMS even for Jlw < Jcrit. Al¬ 
though the concept of one universal threshold Jcrit is ques¬ 
tionable in the formation scenario of SMBH seeds, it is a 
convenient quantification to study the influence of different 
physical processes on the direct collapse scenario. 

4.1 Stellar spectrum 

One should keep in mind that we use the T5 spectrum 
to study the effect of H2 self-shielding (section 12.211 . This 
is an accurate approximation for stars with a character¬ 
istic mass of ~ 100 M0, but subsequent stellar popula¬ 
tions are believed to have softer stel lar spectra and might 
therefore yield a lower value of Jcru (ISugimura et al.l 12014 : 
Agarwal fc Khochfarl l2015l : lAgarwal et al.ll2015^ . However. 
Latif et ^ (l2015fl show that the value of Jcrit only weakly 
depends on the adopted radiation spectra in the range 
2 X 10'‘ K ^ Trad s: 10® K. 

4.2 Resolution 

To test the resolution, we also perform one simulation with a 
better mass resolution, fn our standard approach, we resolve 
the local Jeans mass by 66 Voronoi cells, which correspond to 
a spatial resolution of < O.lpc. For the high-resolution run, 
we increase the mass resolution by a factor of eight, which 
doubles the spatial resolution. The comparison can be seen 
in Fig. 1121 The run with a higher spatial resolution yields 
the same temperature profile as the run with the normal 
resolution. Hence, we can conclude that our normal resolu¬ 
tion is high enough to properly resolve the collapse, and our 
results are not sensitive to the numerical resolution. 

4.3 Photochemistry 

The I Wolcott-Green et all (I2OIII I self-shielding function that 
we use here is intended for use when the H2 is rotationally 
hot and not only the lowest rotational levels are populated. 
This is a reasonable approximation for densities n ~ 10®- 
10 "^cm~®, but at lower densities we would expect most of the 
H2 molecul es to be in the J = 0 o r J — 1 levels, and in this 
regime the IWolcott-Green et al.1 (I2OIII I self-shielding func¬ 
tion will underestimate the effectiveness of H2 self-shielding. 
The effect of this on Jcrit is unclear, and we intend to in¬ 
vestigate this further in future work. In addition, our cur¬ 
rent treatment of the shielding of H2 does not account 
for absorption by the Lyman series lines of atomic hydro - 
gen dHaiinan et al.1 [l9^ . I Wolcott-Green fc HaimanI (l201ll 'l 
show that shielding of H2 by atomic hydrogen becomes im¬ 
portant for column densities of Ah > 10®®cm“®. Including 
this effect may yield larger values for Jcrit and may hence 



log(r/pc) 


Figure 12. Temperature profiles for halo C using the Jeans ap¬ 
proximation, where the green long-dashed curve represents an 
additional run with a two times higher spatial resolution. This 
should be compared to the red short-dashed curve, which repre¬ 
sents a run with the normal resolution and the same background 
flux. The blue dotted curve shows the profile for the same halo 
with normal resolution but a higher LW flux. All curves show 
an isothermal collapse, but the differences induced by a different 
background radiation are larger than those induced by a higher 
resolution. 

exacerbate the difference between the Jeans approximation 
results and the results derived using treecol. 


5 CONCLUSIONS 

We have implemented a new method for the determination 
of H2 column densities in the 3D moving mesh code AREPO 
and used it to study the effect of an improved treatment of 
H2 self-shielding on the ‘direct collapse’ scenario. In a com¬ 
parison to the previously used Jeans approximation, we find 
that the effective column densities are generally smaller with 
our new method and the necessary LW background flux to 
suppress efficient H2 cooling is lower by a factor of about 2. 
More precisely, we find Jcrit — 2000 with our new approach 
compared to Jcrit — 4000 with the Jeans approximation. 
The main reason for this difference is the large directional 
dependence of the self-shielding factor that cannot be cap¬ 
tured with one-dimensional methods. Because the detailed 
morphological and kinematic structure of the cloud matters 
a lot for the determination of the effective column density, 
it is also not possible to find a simple correction factor that 
might reprod uce the results based on tr eecol . 

Following If navoshi fc Tanakal (120141') . the density of pos¬ 
sible direct collapse black hole formation sites scales with 
u-dcbh oc J“;® for Jlw > 10®. Consequently, the factor of 
two, by which Jcrit is lower with our new self-shielding ap¬ 
proach, leads to a number density of direct collapse black 
holes in the early Universe that is about 32 times higher 
than previously expected. Although the number of expected 
direct collapse black holes is significantly higher with our 
new method, the value of Jcrit — 2000 is still too high to 
explain the number density of SMBH at redshift 2 ~ 6 
of nsMBH = 10“® Mpc“® (comoving units) only by the 
isothermal direct collapse scenario. Even under optimistic 
assumptions, Jcrit has to be smaller t han 1000 to explain 
the observed number de nsity (see e.g. IPiikstra et al.ll2014l : 
llnavoshi fc Tanakall201 j l. 
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